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We present a variational ansatz for the ground state of a strongly correlated Bose system. This ansatz goes 
beyond the Jastrow-Feenberg functional form and explicitly enforces coordination shells in the structure of the 
wavefunction. We apply this ansatz to liquid helium-4 with a simple three-variable parametrization of the pair 
functions. The optimized wavefunction is found to give an excellent description of the mid-range correlations 
in the fluid. We also demonstrate the possibility to use this ansatz to study inhomogeneous systems. The phase 
separation and free surface emerge naturally in this wavefunction, even though it is constructed of short-range 
two-body functions and does not contain one-body terms. Because no explicit description of the surface is 
necessary, this provides a powerful description tool for cluster states. 
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I. INTRODUCTION 

The interest in the microscopic nature of the ground state 
of liquid ^He has drawn attention for over half a century and 
has shaped the development of many aspects of the quantum 
many-body theory*. The question continues to be on inter¬ 
est, especially as new correlated bosonic systems are becom¬ 
ing the subject of an experiment, including the cold atomic 
gases^“^. An explicit and numerically efficient expression for 
the many-body wavefunction also has a practical use in com¬ 
puter calculations. A good approximation to the ground state 
reduces the numerical costs and improves the statistical accu¬ 
racy of the true ground state results obtained with the diffusion 
Monte Carlo®’^ as well as the path-integral ground state Monte 
Carlo^“** methods. 

The variational ansatz for liquid ^He has followed 
the path of improving the Jastrow-Feenberg form of the 
wavefunction * ^'' ^, 

i< j i< j<k 

( 1 ) 

where N is the number of atoms, and the A:-body correlation 
factors must have proper symmetry under the exchange of 
particles. Because each successive term in Eq. (1) increases 
the the numerical complexity by an additional factor of N, one 
is in practice limited to two- and three-body terms. Limiting 
Eq. (1) to two-body factors results in the Jastrow function*^’*^. 
In an early work, McMillan*® and Schiff and Verlet*^ used a 
Jastrow function with the two-body function U 2 = —{b/r)^. 
Parameter b was determined variationally. The McMillan 
function captures the most significant features of the system 
caused by the core of the interparticle potential and it con¬ 
tinues to be used successfully as a guiding function for pro¬ 
jector Monte Carlo**'*^’*^. Successive improvements in the 
ground state of helium refined the two- and, later, three- body 
factors in the form (1). Published progress on this topic is 
too numerous to cover in any detail here. Relevant to this 
work, we note the addition of the mid-range correlation^**’^* 
which among other things allowed to replicate the first corre¬ 
lation peak of the parr distribution function g(r); the addition 


of long-range terms in the two-body function U 2 that allows 
to account for the long-wavelength zero-point phonons^^’^^; 
the computation of U 2 based on the maximum overlap with 
the true ground-state^"*’^®; and finally, a detailed optimization 
of the pair factors expanded in terms of the pair scattering 
eigenstates^®’^^ which along with the inclusion of the three- 
body factors allowed to account for nearly all the correlation 
energy. The success of the above works came at the expense 
of the increased complexity and and the number of variational 
parameters that are needed to accurately describe the func¬ 
tions The general functional form of Eq. (1), though, re¬ 
mained unchanged^^. 


The development of the shadow wavefunction (SWE) 
methods^®’®** has to a large degree overtaken the development 
of the wavefunction for liquid helium. The SWE allows to ac¬ 
count for the correlations missed by the Jastrow function, and 
results in an excellent description in terms of both energy and 
structure®*’®^ of "*He. Relevant to this work, we notice that 
SWE can support self-bound states of liquid "*He®®. Shadow 
wavefunction accounts for correlations via integrals on auxil¬ 
iary (shadow) variables. The inclusion of the shadow variables 
may be seen as going beyond the Jastrow-Eeenberg form of 
Eq. (1). However, the integrals on the shadow variables must 
be taken numerically by a Monte Carlo scheme, and in this 
sense shadow wavefunction is not explicit, partially limiting 
its adoption in quantum Monte Carlo. 


We will present a variational ansatz for the ground state 
of liquid "*He which is build upon the Jastrow wavefunction 
but goes beyond the general functional form of Eq. (1). This 
ansatz allows to explicitly control the mid-range structure of 
the liquid and results in a stark improvement of the atomic pair 
distribution already with a three-parameter wavefunction. The 
wavefunction is presented in Section II and the computational 
results are shown in Section III. Section IV presents results 
for inhomogeneous systems, followed by a discussion. 
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II. THE COORDINATED WAVEFUNCTION 

A. Variational ansatz 


Ref. 34 can be seen as a generalized symmetrical form of the 
one-body factor; the coordination part of Eq. (2), however, is 
by the same criterion a full A^-body term. 


Our proposed wavefunction consists of a product of the Jas- 
trow function (limiting Eq. (1) to two-body terms) and of the 
additional term to which we refer as the “coordination term”. 
The wavefunction has the following form, 

Yjciri ,... ,riv) = fj ^y2(|ri - r^'|). (2) 

‘<J ‘ i^i 

The factors y 2 {r) must vanish at large distances. At short dis¬ 
tances, y 2 is expected to raise to a constant. 

The effect of the coordination term in (2) can be seen by in¬ 
spection. Suppose the function y 2 {r) vanishes for distances r 
beyond the mean interparticle distance. In this case, y 2 {rij) 
will have significant value only for the pairs of immediate 
neighbors On the other hand, the number of neighbors 
for each atom is limited by the presence of the repulsive core 
and by the Jastrow part of the wavefunction. Thus the over¬ 
all number of non-vanishing terms y 2 {rij) in the system is, 
roughly speaking, fixed. Under such a restraint, the product 
of sums in the coordination part of Eq. (2) is maximized when 
all sums are equal to each other. That is, the non-vanishing 
values ofy 2 {rij) are distributed equally between the products. 
The wavefunction \j/jc, while constructed only of pairwise 
functions, has a “global” property in that it explicitly demands 
that each atom in the system has an equal expected number of 
immediate neighbors. As we will see, this allows to improve 
the mid-range properties of the system independently of the 
Jastrow factor. 


B. Inspiration and origin 

The inspiration for the coordinated wavefunction \f/jc 
comes from the symmetrized Bose-solid wavefunction pro¬ 
posed by Cazorla et al.^^. This symmetrical solid wave- 
function does an excellent work describing quantum Bose 
solid, both variationally^^ and as a guiding function for im¬ 
portance sampling in quantum Monte Carlo simulations of 
Bose solids^^^®. In fact, one will recognize that Eq. (2) is the 
wavefunction of Cazorla et al., except that the site locations 
of a crystalline structure are here replaced by the positions of 
atoms themselves. 

The solid wavefunction of Ref. 34 forces atoms to be lo¬ 
cated in the vicinity of one of the externally specified lattice 
sites, while at the same time imposing the global restraint by 
favoring single site occupancy. In the liquid, the translational 
symmetry is not broken and there are no preferred positions; 
instead, the atoms in (2) are “localized” around their neigh¬ 
bors. As the overlap of atomic cores is prohibited by the Jas¬ 
trow term, this creates the coordination shells. 

An important distinction between xj/jc of Eq. (2) and the 
symmetrized Nosanow-Jastrow wavefunction of Ref. 34 is in 
the nature of the sum-factors. As discussed in Ref. 35, factors 
that bind atoms to the lattice sites in the solid wavefunction of 


C. Separability 

If the particles of the system are divided into subgroups 
separated by a large distance, the wavefunction ought to re¬ 
duce to the product of the wavefunctions for the individual 
subgroups^*^. While such cluster property is obviously satis¬ 
fied by the Jastrow function, it is less transparent for the co¬ 
ordination term. Suppose all particles are divided into two 
groups, or clusters, A and B. Let the corresponding number 
of particles be Na and Nb, Na-\-Nb = N. The distance be¬ 
tween these clusters is sufficiently large such that the func¬ 
tion y 2 vanishes for any pair of particles from across the two 
groups. 


'ii^AJ &B\y2[\ri-rj\)=0. 

In this case, the coordination sum for any particle in a sub¬ 
group reduces to the sum on that subgroup only, and the coor¬ 
dination term separates. 


V/c(n,---,riv) =nL3'2(lu-r;|) 


= nL3'2(k;-'-rl) X nL3'2(|u-r;|) 


\iCA j^i 


( 


nL3'2(|u-r;|) X 


\ icA 


nL 3 ' 2 (|u-rj|) 

''=B jA‘ 

\ jeB 


= Vk:(u; Wcin'; ,■ ■ ■,) 


i'eA 


i"eB 

= Wc{A)Wc{B). 


Thus a wavefunction for the two clusters reduces to the prod¬ 
uct of the wavefunctions for the individual clusters. 


D. Computational complexity 

The evaluation of the coordinated wavefunction of Eq. (2) 
requires the computation of 0{N^) interparticle distances, and 
the overall computational cost also scales as the second order 
in the number of particles. The scaling holds for the applica¬ 
tion of the Hamiltonian and other relevant operators. To see 
this, we write Eq. (2) as 

VJciru- ■■,rN) = , • • •, '"w), 

i<j i 


with 
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In order to compute all N sums Si, one needs to compute 
N{N — l)/2 values of y 2 {\ri — rj\), so long as the sums are 
stored in memory. This is not a taxing requirement, given one 
must in any case store hN atomic coordinates. Once the sums 
are computed, the computation of the product H; St only re¬ 
quires N operations. 

Similar considerations apply to the computation of the 
Hamiltonian and other relevant expressions. For quantum 
Monte Carlo, one generally needs to compute the contribution 

to the local kinetic energy and the “quantum velocity” 

vector . The relation 

^ = v2logv/+(Vlogv/)^ (3) 

allows us to separate the contributions from the Jastrow and 
the coordination terms. The later is labeled below as y/c- We 
also use a label (■)s_( for the f-th spatial dimension correspond¬ 
ing to particle i; that is, 1 < f < D and s spans from 1 to N. It 
is convenient to define vectors v and u. 


1 2Ci,t 

= n Ej'2(Gt) 

(4) 

1 / \ 1 ; — Xi t 

= 2^F2(G7) „ 

2>i rsi 

(5) 

The quantum velocity is obtained by 


Vs,,log V/c = (« + v).5,,. 

(6) 


The second derivative, summed on the spatial dimension, can 
be written as 


^=1 i^s V 


y'lirsi) + - — -y2{rsi) 
f"si 


X 


X 




(7) 


Notice the cancellation between terms from (7) and (6) sug¬ 
gested by Eq. (3). 

Written in the above form, it is clear that the relevant cal¬ 
culations involve the order of operations with storage re¬ 
quirement of only the first order in N. The calculation may 
proceed as follows. First, one loops through N{N —\)/2 pairs 
of particles and computes the sums S. Then the loop is re¬ 
peated, this time summing the contributions to the vectors v 
and u given by Eqs. (4-5) and the contribution to the second 
derivative given by the first sum on the r.h.s. of Eq. (7). To 
complete the calculation of the kinetic energy, one needs to 
perform N additional operations to compute the second sum 
in the r.h.s. of Eq. (7) and to sum the square of the gradient 
vector according to Eq. (3)"'^*. Thus the computations with the 
coordinated wavefunction of Eq. (2) scales only as the second 
order in the number of particles, although the usual loop over 
the particle pairs needs to be repeated twice. 
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FIG. 1. Variational energy (per particle) for the coordinated wave- 
function given by Eqs. (2), (8)-(10) (bullets) and the Jastrow func¬ 
tion with the McMillan pair factor given by Eq. (8) (triangles). Both 
energies are shown as a function of parameter b which enters the 
two-body correlation factors. For each value of b, the coordinated 
function was optimized with respect to its parameters m and 8. For 
the unit of distance, we use a = 2.556A. The energies were com¬ 
puted for 512-particle systems interacting with Aziz-lE^ pairwise 
potential. Statistical en'ors are smaller than the symbol size. 


E. The form of the pair and coordination factors 

To test the coordinated wavefunction of Eq. (2), we have de¬ 
cided to limit the Jastrow term to the simple McMillan form*® 
with M 2 = —{b/ry. As this term aims to capture the short- 
range correlations in the fluid, the mid-range correlations are 
left to be treated with the coordination term. Having only 
one variational parameter in the Jastrow product simplifies 
the parametrization of the wavefunction. However, the sim¬ 
ple form of the McMillan factor misses over 1 K of the cor¬ 
relation energy, most of it due to its imperfection at short dis¬ 
tances. One should not hope to recover this energy with any 
improvement to the mid-range correlations. 

For calculations, we used the following form of the pairwise 
functions. 


-2{r) = -U-] -T 


yiir) = 1 - exp 


2Lc - r 

m-\ 


R{r) 


Rir) = — 


ir/Lc)' 


4 ’ 



(9) 

( 10 ) 


where b, m, and 5 are the three variational parameters, and 
Lc is the cutoff distance of the calculation given by half the 
dimension of the simulation box. 

The coordination function y 2 was chosen to provide a rea¬ 
sonably sharp cutoff beyond a certain distance 5. At the same 
time, we found it quite important to have a “flat” y 2 at small 
distances, as otherwise the derivative of y 2 interferes with the 
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energy terms produced by the derivatives of the pair factors 
M 2 - Because of this effect , using y 2 of a Gaussian or expo¬ 
nential form results in wide flat energy plateaus in the space 
of variational parameters. Instead, the form given by Eq. (9) 
assures that y 2 reaches a constant at small distances. The rel¬ 
evant small distances are given by the parameter b, and the 
condition can be formulated as 

exp[-(5/Z7)'"] < 1. 

Satisfying the above condition effectively decouples the op¬ 
timization of U 2 and y 2 , allowing for a clear interpretation of 
both terms and for a straight-forward variational optimization. 
Indeed, we found that variationally optimized parameters b, 5 
and m fulfill the above condition to about 10^^. 

As is beneficial for a variational calculation, both Jastrow 
and the coordination terms are symmetrized to result in zero 
gradient of the wavefunction at the computational cutoff Lc- 
The pair term in Eq. (8) is symmetrized in the traditional man¬ 
ner, while the coordination factor y 2 employs a scaling func¬ 
tion R{r) to assure that y 2 vanishes smoothly at the cutoff dis¬ 
tance r = Lc- The use of the scaling function allows for a 
robust implementation of the cutoffat Lc, yet introduces mini¬ 
mal disturbance to y 2 at the relevant distances r Ri 5, as in our 
case < 10^^. We found that using scaling function 

provides a convenient way for symmetrizing the wavefunc¬ 
tion. 


III. RESULTS 
A. Variational optimization 

We carried the variational optimization with three- 
dimensional 512-atom ^He system at the equilibrium 
densityof liquid ^He, po = 0.365(7^^ = 21.8 nm^^. Here 
and below, we use the reduced unit of length equal to O’ = 
2.556 A. All observables where computed on Markov chains 
generated by the Metropolis method^"*^ with single-particle up¬ 
dates. We used a GPU cluster to speed up the calculations us¬ 
ing a modification of the QL quantum Monte Carlo package^^. 


TABLE I. Optimized parameters of the coordinated and non- 
coordinated Jastrow wavefunctions with the McMillan factor. Dis¬ 
tances are specified in units of o = 2.556A. Lowest line shows the 
thermodynamic limit extrapolation of the per-particle energy, with 
up to 1920 particles used for the calculation. The interaction was 
modeled with the Aziz pair potential from Ref. 42. 



Coordinated 

Non-coord. 
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— 
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FIG. 2. Pair distribution function g[r) as a function of the interparti¬ 
cle distance, obtained for a 512-atom system. For the unit of distance, 
we use a = 2.556 A. Unconnected black errors bars: unbiased (pure) 
estimator obtained with DMC, as described in the text. Connected 
green triangles: Jastrow function with the McMillan pair factor as 
specified in Eq. (8). Connected red bullets: energy-optimized three- 
parameter coordinated wavefunction \jfjc with the McMillan factor, 
given by Eqs. (2),(8)-(10). Errors bars for both VMC calculations 
are smaller than their corresponding symbol sizes. The inset shows 
the details of the first correlation peak. 

The coordinated wavefunction of Eq. (2) was taken in three- 
parameters form given by Eqs. (8-10). The system Hamilto¬ 
nian 

' KJ 

was used with the pair potential by Aziz^^. 

The three variational parameters b, m, and 5 were opti¬ 
mized on a grid. Eigure 1 shows variational energy as a func¬ 
tion of parameter b, given optimal m and 5 for each value of 
b. The results are compared to the (non-coordinated) Jastrow 
function with the McMillan factor given by Eq. (8). Optimal 
value of the parameter b for for the coordinated function was 
found toheb= 1.19(7, slightly below the optimal value of the 
non-coordinated function, b = 1.20(7. As expected, we found 
little variation in optimal value of 5 with respect to chang¬ 
ing the value of parameter b. In the range shown in Eigure 1, 
optimal 5 varies less than two percent. We also notice rela¬ 
tively weak correlation between parameters 5 and m near the 
variational minimum. 

The optimized values of variational parameters are shown 
in Table I. The table also shows the value of optimized energy 
extrapolated to the thermodynamic limit. Eor comparison. Ta¬ 
ble I also lists the energy for the Jastrow function with the 
McMillan factor. This energy differs slightly from the one ob¬ 
tained by McMillan*®, which can be prescribed to the differ¬ 
ence in the interaction potential. As expected, one will notice 
that the gain in the correlation energy is mild, and amounts 
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to just under 200 mK. This is in part due to the fact that the 
missing mid-range correlations are not responsible for a large 
amount of energy, but also because the presence of the coordi¬ 
nation term ever so slightly offsets the correlation hole which 
in turn carries an energy penalty. 


B. Structural properties of the coordinated function 

As both the potential energy and the wavefunction are built 
from the pairwise functions, the properties of the system are 
captured by the pair distribution function. The computed 
pair distribution function g{r) is shown in Figure 2, along 
with the results for the McMillan function and an unbiased 
(pure) estimate for the g{r) obtained with the diffusion Monte 
Carlo (DMC). The unbiased DMC estimator for g{r) was ob¬ 
tained with the ancestry tracking algorithm of Casulleras and 
Boronat^®. Such an unbiased estimator is computed from the 
projected ground state and can be expected to reflect accu¬ 
rately on the experimental values'^®”'^^. In properly converged 
calculations, pure DMC results do not depend on the DMC 
guiding function. However, it is worth pointing out that the 
guiding function for the DMC calculation was in fact the Jas- 
trow function with the McMillan factor and it did not contain 
the coordination factor. In all three cases shown in Figure 2, 
the calculations were performed with 512-particle systems at 
the equilibrium density of ^He, p = 0.365(7 = 21.8 nm^^. 
The variational parameters were chosen by energy optimiza¬ 
tion, as specified above, and are given in Table I. 

It is notable that the coordinated wavefunction reproduces 
accurately the first correlation peak in the pair distribution 
function. The inset in Figure 2 shows the detail of the 
first maximum. The first correlation minimum is reproduced 
slightly less accurately. The following oscillations in the pair 
distribution function are also reproduced better by the coordi¬ 
nated wavefunction, albeit with decreasing accuracy. The po¬ 
sition of the maxima and minima in the pair distribution was 
also considerably improved by the coordination term. The de¬ 
tails are given in Table II. However, the absolute value of 
these successive oscillations is minute, and they are at dis¬ 
tances where the pair potential is vanishing rapidly. Thus their 
influence to the overall energy is nonsignificant. 


TABLE II. The degree to which the computed pair distribution func¬ 
tions g{r) capture the unbiased estimate g*(r). The values are com¬ 
puted as \g(r,„) ~g*(rm)|/|l ~g*(rm)|, where r„, are the locations 
of extrema of g*(r). The simulation conditions are described in Fig¬ 
ure 2. 
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Coordinated 
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0.9 

0.7 

0.6 

0.4 

0.4 

Non-coor. 

0.8 

0.6 

0.5 

0.3 

0.2 

0.2 


IV. INHOMOGENEOUS SYSTEMS 

Jastrow wavefunctions based on a short-range pair factors 
cannot support the formation of a self-bound state. That is, 
a simulation in a sufficiently large box will result in a low- 
density uniform gas with near-zero potential and kinetic en¬ 
ergy. Helium liquid, however, is self-bound. To describe in¬ 
homogeneous systems, one generally adds one-body factors 
which bind the liquid phase to a desired shape^^“^*. This 
has obvious disadvantages if the surface shape is complex, 
and may pose additional challenges when one needs to main¬ 
tain the translational symmetry in the system^^. Parametriza- 
tion of the surface adds to the required number of the varia¬ 
tional variables. Self-binding may also be enforced through 
the use of long-range terms in the two-body factors, such 
as introduced in Ref. 53, with additional term in two-body 
function U 2 {r) proportional to the distance between the parti¬ 
cles r. Such a wavefunction serves well as a trial wavefunc¬ 
tion for a projector Monte Carlo calculation, yet variation- 
ally, kinetic per-particle energy of a system with 1/2 —OCr 
for large r is divergent with the increasing number of particles 
A as ^ [h^/m)py^N^l^a, where po is the bulk density. This 
presents a number of challenges, as at the very least a must 
be A-dependent. 

The coordinated wavefunction has an unexpected feature 
in that by design it supports a self-bound state of the atoms. 
Upon inspection, one will notice that the coordination sum- 
factors in Eq. (2) in fact vanish in the limit of low-density, 
uniformly distributed gas. Thus the coordination term requires 
that atoms form clusters, so far as function y 2 {r) falls off suf¬ 
ficiently rapidly with distance. The size and number of the 
clusters is determined by the variational parameters and parti¬ 
cle density. For instance, a gas of dimers already has non-zero 
coordination term. Increasing the range of y 2 (which in our 
case translates to increasing 5 or decreasing m) increases the 
size of the clusters. The parameters also provide control over 
the structure of the surface. 

We have carried variational calculation with the coordi¬ 
nated wavefunction with 1000 atoms in a simulation box that 
resulted in particle density 10^^ (7^^ w 0.06 nm^^. The Jas¬ 
trow wavefunction results in a ground state of dilute gas with 
close to zero energy. However, the coordinated wavefunction 
resulted in a bound state for a wide range of parameters m and 
5. Only small 5 resulted in the unbound states albeit with 
positive energy. We also find that the extent of function y 2 
controls the average cluster size and thus the energy. Varia¬ 
tional optimization of m and 5 results in a state with a single 
liquid droplet. 

To demonstrate the robustness of the inhomogeneous sim¬ 
ulation, we carried Monte Carlo sampling of the coordinated 
wavefunction with parameters m = 6.55 and 5 — 4.50(7 (i.e., 
with y 2 having larger extent than for the bulk). The initial 
coordinates of the 1000 particles were randomly distributed 
in the simulation box. The sampling sequence is presented 
in Figure 3. Soon after the start, the Markov chain arrives at 
configurations with multiple small clusters. As the clusters 
merge, a single droplet is eventually formed. The center of 
mass of the system is not fixed, and the droplet continues to 
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FIG. 3. Progression of Markov chain during Metropolis sampling of a system with N = 1000 atoms with the coordinated wavefunction. 
The size of the cubic box is equal to L = 100 a = 256 A, which would correspond to a very dilute homogeneous system. Top row, left to 
right, shows the initial state of the system (with atoms distributed randomly and uniformly), and the system correspondingly after 10^, 10^, 
and 10^ macroupdates. After 10^ updates, the Markov chain reaches droplet configuration which is then sampled throughout the (periodic) 
simulation cell. The bottom row, left to right, shows configurations after 2- 10^, 3 • 10^, 10^, and 2- 10^ macroupdates. Wavefunction 
parameters b = 1.19a, S = 4.60a, m = 6.55. Metropolis sampling was carried via single-particle updates, with a fixed Gaussian-distribution 
of displacements which resulted in the acceptance ratio of above 20% in the homogeneous phase to below 35% in the condensed phase. Each 
“macroupdate” equals N single-particle Metropolis attempts. 


sample the entire simulation cell. 

The inner structure of the droplets and clusters depend 
strongly on the two-body function 112 - However, Jastrow func¬ 
tion with the McMillan factor underestimates the equilibrium 
density of the bulk ^He. Without the fixed density constraint, 
it is to be expected that the inhomogeneous simulation should 
result in lower densities of the condensed phase. This was in¬ 
deed observed. For example, the droplet shown in Figure 3 
has inner density that is less than 70% of the bulk equilibrium 
helium density. Thus the droplet calculation presented here 
should be seen as a demonstration of principle. The details 
of their structure, which require a more detailed Jastrow term, 
will be the subject of further investigation. 


V. CONCLUSION 

We have considered a wavefunction ansatz for strongly cor¬ 
related Bose system that goes beyond the Jastrow-Feenberg 
expansion. Originating from a symmetrical solid wavefunc¬ 
tion proposed by Cazorla et al.^'^, it is a Bose-liquid wave- 
function which explicitly promotes the creation of the coor¬ 
dination shells around atoms. The function is translationally 
and exchange symmetric. It is fully explicit and is computa¬ 
tionally hard as 0{N^), making it well suitable for treatment 


with quantum Monte Carlo. 

To demonstrate the coordination effect, we have studied the 
wavefunction with the one-parameter McMillan factor for the 
Jastrow term, and a two-parameter coordination function. The 
resulting three-parameter wavefunction was straight-forward 
to optimize variationally. The short-range nature of the 
McMillan factor allowed to directly observe the effects of 
the coordination terms on the mid-range structure of the liq¬ 
uid. Indeed, the optimized wavefunction results in superior 
description of mid-range correlations in the system. Compar¬ 
ing with unbiased estimate for the pair distribution function 
obtained with the diffusion Monte Carlo, we find that the first 
correlation peak is reproduced almost exactly. Moreover, the 
structure of the pair distribution function is improved consis¬ 
tently throughout larger distances as well. 

As was first demonstrated in Ref. 20, the first correlation 
peak can be reproduced rather exactly with the Jastrow func¬ 
tion. However, this required eight variational parameters, and 
already the description of the first minimum was significantly 
lacking. Other approaches to accurately describe the mid¬ 
range structure with the Jastrow factors alone have also been 
reported^^. In our case, the addition of the coordination term 
allows to separate the short- and middle-range correlations, 
which can be accounted for correspondingly by the Jastrow 
and the coordination terms. 
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By construction, the coordinated wavefunction supports a 
self-bound state. Consequently, the simulation of inhomo¬ 
geneous systems does not require the addition of one-body 
terms. Moreover, inhomogenuity and surface formation at low 
densities result directly from the variational optimization of 
the bulk wavefunction. Since the variational ansatz does not 
require knowledge of the surface geometry, this also provides 
a powerful tool for cluster states of matter. However, we find 
that a satisfactory description of the inhomogeneous phase of 


helium requires improvements in the Jastrow pair term, which 
was here limited to the McMillan form for simplicity. 

The separation of the mid-range correlations into the co¬ 
ordination term which was demonstrated here means that the 
Jastrow pair term in the coordinated wavefunction only needs 
to account for the short-range correlations and possibly for the 
well-understood long-range correlations arising from zero- 
point phonons. This makes it promising that an accurate 
short-range pair term can be designed in future with a sim¬ 
ple parametrization. 
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